#include "LinerEqus.h"

LinerEqus::LinerEqus()
{
    row = col = 0;
    m_a = m_b = m_x = l = NULL;
    true_flag = 0;
}

LinerEqus::LinerEqus(int M, int N)
{
    row = M;
    col = N;
    //m_ivec.clear();
    m_a = SetMatrix(M, N, 0);
    m_b = SetMatrix(M, 1, 0);
    m_x = SetMatrix(M, 1, 0);
    l = SetMatrix(M, N, 0);

    true_flag = 0;

}

LinerEqus::~LinerEqus()
{
    FreeMatrix(m_a, row, col);
    FreeMatrix(m_b, row, 1);
    FreeMatrix(m_x, col, 1);
    FreeMatrix(l, row, col);
}

//动态建立M*N的数组，并输入值。返回建立的数组的地址
double** LinerEqus::SetMatrix(int M, int N, bool input_flag )
{
    double **p = NULL;
    p = (double **)malloc(M * sizeof(double *) );
    for (int i = 0; i < M; i++ )
    {
        p[i] = (double *)malloc(N * sizeof(double) );
    }
    if(!input_flag)
        return p;                                                   //如果input_flag为false，返回p
    cout << "Please enter the value of the elements:" << endl;      //手工输入矩阵元素的值
    for(int i = 0; i < M; i++)
    {
        //cout << "请输入第" << i+1 << "行数据：" << endl;
        for(int j = 0; j < N; j++ )
        {
            cin >> p[i][j];
            //cout << p[i][j] << endl;
        }
    }
    return p;
}

//打印矩阵
void LinerEqus::PrintMatrix(double **p, int M, int N)
{
    for(int i = 0; i < M; i++ )
    {
        for(int j = 0; j < N; j++ )
            cout << p[i][j] << '\t';
        cout <<endl;
    }
}

//释放动态数组p
void LinerEqus::FreeMatrix(double **p, int M, int N)
{
    for(int i = 0; i < M; i++ )
    {
        free(p[i]);
        p[i] = NULL;
    }
    free(p);
    p = NULL;
}

void LinerEqus::LUDec(double **a, double **b, int M, int N)
{
    int t;
    for(int i = 0; i < M; i++)
    {
        t = FindMaxCol(a, M, N, i);
        m_ivec.push_back(t);					//将最大主元所在的行保存下来；
        if (t != i )
        {
            ExcRow(a, M, N, i, t);
            ExcRow(b, M, 1, i, t);
        }
        RowOp(a, b, M, N, i ); //消元
    }
}

//找出数组p的第k+1列中，在第k+1行（包括）后绝对值最大的元素,返回绝对值最大的元素所在的行（矩阵中第i+1行）
int LinerEqus::FindMaxCol(double **p, int M, int N, int k)
{
    double temp = p[k][k];
    int i = k;
    for(int j = k; j < N; j++)
    {
        if(fabs(temp) < fabs(p[j][k]) )
        {
            i =j;
            temp = p[j][k];
        }
    }
    return i;
}

//交换矩阵p的第k+1行和第t+1行
void LinerEqus::ExcRow(double **p, int M, int N, int k, int t )
{
    double temp;
    for(int j = 0; j < N; j++)
    {
        temp = p[k][j];
        p[k][j] = p[t][j];
        p[t][j] = temp;
    }
}

//初等行变换
void LinerEqus::RowOp(double **a, double **b, int M, int N, int k)
{
    double temp;
    for(int i = k + 1; i < M; i++)
    {
        temp = a[i][k] / a[k][k];

        if(true_flag)
            l[i][k] = temp;

        b[i][0] -= temp * b[k][0];
        for(int j = 0; j < N; j++ )
        {
            a[i][j] -= temp * a[k][j];
        }
    }
}

//迭代求出解
void LinerEqus::Solve(double **a, double **b, double **x, int M, int N)
{
    double sum;
    x[M-1][0] = b[M-1][0] / a[M-1][M-1];
    for(int i = M - 2; i >= 0; i--)
    {
        sum =0;
        for(int j = i + 1; j < M; j++)
            sum += a[i][j] * x[j][0];
        x[i][0] = (b[i][0] - sum ) / a[i][i];
    }
}

//将数组sou的数据拷贝至数组des
void LinerEqus::CpArray(double **sou, double **des, int M ,int N )
{
    for(int i = 0; i < M; i++ )
    {
        for(int j = 0; j < N ; j++ )
        {
            des[i][j] = sou[i][j];
        }
    }

}


//根据m_ivec中值的情况，先对a、b进 行变换 ，然后进行LU分解
void LinerEqus::LUTure(double **a, double **b, double **x, int M, int N )
{
    vector<int >::iterator iter = m_ivec.begin();
    for(int i = 0;iter != m_ivec.end(); iter++ )
    {
        if (*iter != i)
        {
            ExcRow(a, M, N, i, *iter);
            ExcRow(b, M, 1, i, *iter);
        }
        i++;
    }
    true_flag = 1;
    //初始化L矩阵
//    if(l == NULL )
//        l = SetMatrix(M, N, 0);
    for (int i = 0; i < M; i++)
        for (int j = 0; j < N; j++)
            if (i == j )
                l[i][j] = 1;
            else
                l[i][j] = 0;


    LUDec(m_a, m_b, M, N);
    Solve(m_a, m_b, m_x, M, N);

}
